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SUMMARY 


System dynamics approach offers great potential for addressing how intervention policies can 
affect the spread of emerging infectious diseases in complex and highly networked systems. Here, 
we develop a model that explains the severe acute respiratory syndrome coronavirus (SARS-CoV) 
epidemic that occurred in Hong Kong in 2003. The dynamic model developed with system 
dynamics methodology included 23 variables (five states, four flows, eight auxiliary variables, 

six parameters), five differential equations and 12 algebraic equations. The parameters were 
optimized following an iterative process of simulation to fit the real data from the epidemics. 
Univariate and multivariate sensitivity analyses were performed to determine the reliability of 
the model. In addition, we discuss how further testing using this model can inform community 
interventions to reduce the risk in current and future outbreaks, such as the recently Middle 
East respiratory syndrome coronavirus (MERS-CoV) epidemic. 
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INTRODUCTION with the severe acute respiratory syndrome 
(SARS)-CoV that produced a global epidemic with 
more than 8000 human cases in 2002-2003 [8, 9]. 
First, a number of patients infected with both viruses 
developed an acute respiratory disease that in some 
cases resulted in death [7, 9]. In this sense, 
MERS-CoV appears to be highly pathogenic with 
an estimated case-fatality rate of around 50%, 
although this might be an overestimation as many 
infected patients may not have sought hospital assist- 
ance [2]. Second, both MERS-CoV and SARS-CoV, 
* Author for correspondence: Dr E. Alvarez, Centro de Biologia belong to the genus Betacoronavirus and are closely re- 
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Middle East respiratory syndrome (MERS), is a res- 
piratory illness caused by a novel coronavirus (CoV) 
[1]. The disease was reported for the first time in 
Saudi Arabia in June 2012 and spread to several coun- 
tries in Africa, Asia, Americas and Europe [2, 3]. The 
capability of human-to-human transmission has been 
observed in at least four hospital settings [4-7]. 
Significantly, MERS-CoV shares certain similarities 
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through intermediate species (e.g. camels for 
MERS-CoV and civet cats for SARS-CoV). Third, 
since they are new emerging viruses, there are no effec- 
tive vaccines or antiviral treatments. The SARS epi- 
demic is a clear example of how a networked health 
system can respond to a new threat to human health. 

One of the best-characterized outbreaks during the 
SARS epidemic was in Hong Kong in 2003 where 
there were 1755 confirmed cases with 299 deaths 
(WHO; http://www.who.int/csr/sars/country/en/index. 
html). The outbreak began in mid-February caused 
by an infected person who travelled from 
Guangdong to Hong Kong [13]. An important fact 
in the generation of a model of the outbreak is that 
the Hong Kong health authorities quickly implemen- 
ted contagion control procedures [14]. In general, two 
interventions were introduced to prevent the spread of 
SARS-CoV. The first was implementation of quaran- 
tine measures to isolate healthy people who had been 
in contact with infected people and therefore poten- 
tially in contact with the virus; isolating those that 
could be infected and asymptomatic during the incu- 
bation period and isolating and treating patients 
who had developed the disease. The other intervention 
was the application of protective measures by healthy 
people who were in contact with infected people to 
avoid becoming infected, such as respiratory protec- 
tion for healthcare workers and daily disinfection of 
the environment of affected rooms [14]. These control 
interventions were implemented progressively in the 
Hong Kong Special Administrative Region from 
mid-March to late April [15]. The application of 
these procedures allowed the rapid control of the out- 
break in the subsequent months. In this sense, it is esti- 
mated that the epidemic in Hong Kong ended in late 
June. 

System dynamics has been proved to be a powerful 
instrument for analysing social, economic, ecological 
and biological systems [16]. In addition, disease epi- 
demiology has been studied using this approach, 
whereby system dynamics offers the practical appli- 
cation of concepts by computerized models that 
allow the systematic test of different scenarios and 
alternative policies [17-19]. In this work, we per- 
formed a modelling of the Hong Kong SARS-CoV 
outbreak using system dynamics. The developed 
model contains five states, four flows, eight auxiliary 
variables and six parameters that interact through 
five differential and 12 algebraic equations. The para- 
meters of the model were optimized following an iter- 
ative process of simulation to obtain a model that 


largely fits the data available to the epidemic. 
Moreover, the credibility of our model and its par- 
ameters are supported by both univariate and multi- 
variate sensitivity analyses. The model reproduces 
how the implementation of control measures was ef- 
fective in preventing the spread of infection to the 
rest of the population. Basically, these measures result 
in a sustained reduction in the frequency of contacts. 
At present, the application of similar measures for in- 
fection containment can help to prevent the spread of 
new emerging epidemics, such as the outbreak caused 
by MERS-CoV. 


METHODS 
Data sources and collection 


Data on the total population of Hong Kong in 2003 
was obtained from statistics of the Census and 
Statistic Department, The Government of the Hong 
Kong Special Administrative Region (http://www.cen- 
statd.gov.hk/home/). Data on cumulative cases, 
deaths and recoveries during the SARS epidemics in 
Hong Kong was obtained from Global Alert and 
Response databases of the WHO (http://www.who. 
int/csr/sars/country/en/index.html). 


Modelling 


The model was developed following the four-step se- 
quence proposed by system dynamics methodology 
[16]. First, the real data from the Hong Kong SARS 
epidemics (Fig. 1) together with other evidences and 
our professional experience were used to create a 
mental modelling of the reality of the outbreak. 
Second, the model structure that is able to explain 
the evolution of the epidemics was represented as a 
Forrester diagram (Fig. 2). Third, the outbreak was 
mathematically modelled as a continuous dynamic 
process represented by a set of differential and al- 
gebraic equations (Tables 1-3). Finally, the model 
was optimized to fit the real data from Figure 1. 

A dynamic compartmental model provides a frame- 
work for the study of transport between different com- 
partments of a system, i.e. well known epidemiological 
compartmental models [20]. To explain how the pro- 
tective measures taken by the government of Hong 
Kong allowed the rapid control of the epidemic we 
consider a dynamic model with five compartments 
(states) and four transitions (physical flows) between 
them. This model is based on two assumptions. 
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Fig. 1. Real data from the Hong Kong SARS outbreak. (a) Daily reported new cases, deaths, and recovered. (b) Cumulative 
cases, dead, and recovered. 
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Fig. 2. Graphical description of the model using Vensim software. Stock variables are represented inside boxes; flow 
variables are designed as arrows with a valve; the rest of the variables are auxiliary (inside circles) or parameters (in bold); 
single lines with arrow are used as connectors to specify that the destination variable is affected by the variable of origin. 


every day there is a different number of people: who 
are infected without symptoms of the disease (inci- 
dence); who develop signs and symptoms of the illness 
(sick per day); who recover from the disease 


First, the individuals are classified into five subgroups 
(SUSCEPTIBLE, LATENT, INFECTED, RECOVERED, DEAD). 
Although, the last subgroup is not strictly needed, it 
is used to keep an account of those dead. Second, 
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Table 1. Differential equations for the five model stocks 


Eqn no. Stock Equation Initial value Units 

(1) SUSCEPTIBLE 0 SUSCEPTIBLE(f)/0t= —Incidence(t) 6:5 x 10°* People 

(2) LATENT 0 LATENT(‘))/ot = Incidence(¢) — Sick per day(?) 0 People 

(3) INFECTED 0 INFECTED(t)/ot = Sick per day(t) — Daily lf People 
recovered(t) — Daily deaths(?) 

(4) RECOVERED 0 RECOVERED(t)/Ot = Daily recovered(t) 0 People 

(5) DEAD 0 DEAD(t)/dt = Daily deaths(t) 0 People 


* Hong Kong’s population in 2003 was about 6:5 million people who were susceptible to infection with the SARS virus. 
+ It is considered that the Hong Kong outbreak originated from an infected person from the province of Guangdong [13]. 


Table 2. Algebraic equations for the four model flows 


Eqn no. Flow Equation Units 

(1) Incidence Incidence(t) = SUSCEPTIBLE(t) X Prevalence(t) x Contagion rate(t) People day! 

(2) Sick per day Sick per day(t) = LATENT(‘)/Incubation period People day~' 

(3) Daily recovered Daily recovered(t) = INFECTED(f) X Recovery rate(t) People day~! 

(4) Daily deaths Daily deaths(t) = INFECTED(‘) X Case fatality/Disease duration People day~! 

Table 3. Algebraic equations for the eight model auxiliary variables 

Eqn 

no. Auxiliary variable Equation Units 

(1) Prevalence Prevalence(t) = INFECTED(f)/Population(t) Dimensionless 

(2) Recovery rate Recovery rate(¢) = (1 — Case fatality)/(Disease duration) Day~' 

(3) Contagion rate Contagion rate(t) = Frequency of contacts(t) x Infectivity Day! 

(4) Population Population(t) = INFECTED(t) + LATENT(f) + RECOVERED(f) + SUSCEPTIBLE(f) People 

(5) Cumulative cases Cumulative cases(t) = DEAD(t) + RECOVERED(f) + INFECTED(f) People 

(6) Cumulative attack rate | Cumulative attack rate (t) = Cumulative cases(f)/Population(s) Dimensionless 

(7) Frequency of contacts Frequency of contacts(t) = Daily contacts /(1 + (Cumulative attack rate/_ Day~! 
Threshold of cumulative attack rate)*) 

(8) Basic reproductive Basic reproductive number(t) = Contagion rate(t)/Recovery rate(¢) Dimensionless 


number (Ro) 


(daily recovered), and who die (daily deaths) as conse- 
quence of the outbreak. 

The model structure built with Vensim software 
(Ventana Systems Inc., USA) (Fig. 2) contains the 
five states and the four flows mentioned above and 
also eight auxiliary variables (inside circles) and six 
parameters (in bold). These elements are linked by 
the physical flows (double line with arrow) and by 
the information transmissions (single line with 
arrow) according to the mathematical model repre- 
sented by the set of differential and algebraic equa- 
tions (Tables 1-3). 

The five differential equations in Table | establish 
the mass balance (inflows minus outflows) in the 


compartments, and as such, they describe the changes 
in the number of people (stocks) in the five subgroups. 
The four algebraic equations of Table 2 express that 
the physical transitions depend directly from the 
stock in the compartment of origin and indirectly 
from other stocks and parameters through the corre- 
sponding auxiliary variables. These equations involve 
three stocks (SUSCEPTIBLE, LATENT, INFECTED), three 
auxiliary variables (prevalence, contagion rate, recov- 
ery rate) and three parameters (incubation period, 
case fatality, disease duration). For instance, the 
incidence flow is proportional to the number of sus- 
ceptible people, the prevalence and the contagion 
rate. The contagion rate or transmission coefficient 


theoretically depends on the number of contacts per 
unit time and the probability of effective contact, 1.e. 
the probability that a contact between an infectious 
source and susceptible host results in a successful 
transfer whereby the susceptible host becomes infected. 

The first three algebraic equations of Table 3 ex- 
press the auxiliary variables (prevalence, contagion 
rate, recovery rate) used in the equations of Table 2. 
In turn, these depend on other auxiliary variables, 
stocks and parameters. For instance, prevalence is 
defined as the ratio between number of infected people 
and the total population [Table 3, equation (1)]. The 
recovery rate depends on the illness duration and 
case fatality [Table 3, equation (2)]. Equation (3) in 
Table 3 expresses the contagion rate which is directly 
dependent on the auxiliary variable ‘frequency of con- 
tacts’ and the parameter ‘infectivity’. 

One of the key variables in the model is frequency 
of contacts because it tries to reproduce the control 
measures carried out by the Hong Kong Government 
to control the SARS outbreak. We assumed that: (1) 
the control measures gave, as consequence, a marked 
reduction of daily contacts in Hong Kong; (2) the con- 
trol measures were based on the cumulative attack 
rate, measured as the ratio between cumulative cases 
and total population [Table 3, equation (6)]. There- 
fore, we decided to use the frequency of contacts as 
the daily contacts modulated by the repression Hill 
function (Fig. 3), in accordance with the equation 
(7) of Table 3. Despite the complexity level of biologi- 
cal systems several cases have been modelled using the 
Hill function, in order to simulate repressor activities 
of enzymatic reactions and the regulation mechanisms 
of several transcription factors [21]. In our model this 
repressor function allows the establishment of the 
relationship between frequency of contacts and the 
cumulative attack rate. Note that the cumulative 
attack rate is used as an active repressor, so half- 
maximal repression occurs when the cumulative 
attack rate is equal to the threshold, and almost 
total repression occurs when this cumulative attack 
rate is double the threshold. 

Equation (8) in Table 3 is used to report on the 
basic reproductive number (Ro), which is defined as 
the expected number of secondary infectious cases 
generated by an average infectious case in an entirely 
susceptible population [20]. In our model this number 
is calculated as the ratio between the contagion rate 
and the recovery rate. If Ro < 1, then the infected indi- 
viduals in the total population fail to replace them- 
selves, and the disease does not spread. However, if 
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Fig. 3. Normalized representation of repression Hill 
function. The solid line corresponds to a value of n =3 [as 
used in equation (7) of Table 3]. The dotted line 
corresponds to the threshold repression function (when n 
tends to infinity). The representation is normalized to the 
threshold of cumulative attack rate on the abscissa axis 
and the value of the discontinuity (daily contacts) on the 
ordinate axis. 


Table 4. Values for the six model parameters 


Parameter Value Units 
Infectivity 0-022 Dimensionless 
Incubation period 533 Day 

Disease duration 26 Day 

Case fatality 0-17 Dimensionless 
Daily contacts 16-8 Day"! 
Threshold of cumulative 7-8 x 10> Dimensionless 


attack rate 


Ro>1, the number of cases generally increases over 
time and the disease spreads. 


Parameterization 


One of the most challenging issues in system dynamics 
modelling is to establish the value of the model par- 
ameters. The parameters can be estimated taking ad- 
vantage of the known information available in the 
literature and can be optimized by means of an itera- 
tive process of simulation. Using this approach, we set 
values for the six parameters of our model which are 
summarized in Table 4. 

Infectivity expresses the ability of the pathogen to 
penetrate, survive and multiply in the host and it is 
measured through the secondary attack rate which is 
defined as the probability that infection occurs in sus- 
ceptible persons within a reasonable incubation period 
following known contact with an infectious person or 
an infectious source. Epidemiological studies in 
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Table 5. Probabilistic distribution patterns of the model parameters 


Parameter Distribution Range Distr. values* Reference 
Infectivity Gamma (0-001-0-2) 1-1/0/0-025 [23] 
Incubation period Gamma (1-30) 1-4/1/3-85 [24] 
Disease duration Gamma (1-70) 5-6/0/4:55 [24] 

Case fatality Normal (0-14-0-21) 0-171/0-01 [25] 
Daily contacts Normal (1-40) 16-8/8-5 [29] 
Threshold of Normal (5-0 x 10->-1-1 x 1074) 7-8 x 10-91 x 10-> 


cumulative attack rate 


* The distribution values for normal distributions are mean and standard deviation, whereas for gamma distributions they are 


order/shift/stretch. 


Singapore showed that 80% of individuals infected 
with SARS-CoV did not cause secondary infections, 
suggesting low infectivity [22]. After the optimization 
process, the final value set for infectivity was 0-022, 
which is in agreement with a previous work examining 
the probability of transmission for SARS-CoV [23]. 

The incubation period of the disease, during which 
time the individual is asymptomatic appears to be be- 
tween 2 and 10 days with an average value of 5-3 days 
[24]. Similarly, although the duration of the disease is 
variable depending on the case, the average of the 
most severe and the mildest cases is 26 days [24]. 
During this time period the virus can be transmitted 
to other people. 

Several studies have shown that the case fatality of 
SARS-CoV was also variable in the 2003 outbreak 
depending on different factors. It has been estimated 
that the mean case fatality for SARS in Hong Kong 
was around 0-17 [25]. 

Social contact patterns have been shown to be 
highly relevant on the transmission dynamics of res- 
piratory infections such as measles, rubella and 
influenza [26-28]. Quantifying this parameter is criti- 
cal for estimating the impact of such infections, for de- 
signing and targeting preventive interventions, and for 
modelling their impact. In our model the daily con- 
tacts parameter, indicating the mean number of con- 
tacts of a person in | day, was set to 16-8 after the 
optimization process, which is in agreement with pre- 
vious studies quantifying these social mixing patterns 
[29, 30]. 

The parameter ‘threshold of cumulative attack rate’ 
is critical in our model since it allows setting the value 
of the cumulative attack rate at which control meas- 
ures were established by the Hong Kong authorities. 
This parameter was finally set at 7:8 10~°. 
Therefore, according to the Hill function of 
Figure 3, for any number of the cumulative attack 


rate below the 78 cases by million, the frequency of 
contacts will be greater than 8-4 (half of the daily con- 
tacts), and in the opposite case the frequency of con- 
tacts will be markedly reduced. 

As explained before, our dynamic model was sub- 
jected to successive rounds of simulation and optimi- 
zation in order to fine-tune the parameters of the 
model. In all these simulations we assumed that the 
Hong Kong epidemic originated from a traveller 
from Guangdong in China [13], and at the beginning 
of the outbreak the entire population of Hong Kong 
was susceptible to SARS-CoV infection. These 
assumptions are represented by the initial values in 
the five stocks (Table 1). Moreover, based on the 
reported data, the simulation period was set to 146 
days with a time step of 1 day. 


Sensitivity analysis 


Parameters of system dynamics models are subject to 
uncertainty. Sensitivity analyses were conducted to pro- 
vide insight into how uncertainty in the parameters 
affects the model outputs and which parameters tend 
to drive these variations. In this task it is essential to 
define the probabilistic distribution patterns of the 
model parameters, which are shown in Table 5, together 
with the references that support these patterns. 

The most influential parameters were estimated by 
univariate analyses, in which changes in the model 
output were studied after disturbance in each par- 
ameter value independently. In complex models, uni- 
variate sensitivity analysis can be insufficient for a 
comprehensive study of the model. Simultaneous 
fluctuations in the value of more than one parameter 
may create an unexpected output change due to non- 
linear relationships in different model components. 
Thus, to test the influence of simultaneous changes 
in the model parameters, the univariate analyses 
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Fig. 4. Simulation of the model. For simulation, the time period was set to 146 days with a time step of 1 day. The 
graphs represent the output of the simulation (solid lines) and the data reported for the SARS outbreak (dotted lines). 
(a) Sick per day, (b) cumulative cases, (c) recovered, (d) dead, (e) daily recovered, and (f) daily deaths. 


were followed by Monte-Carlo multivariate sensitivity 
analysis, in which the values of the six parameters 
were altered at the same time. 


Software 


Modelling process, simulations and sensitivity analy- 
ses were performed using Vensim DSS software 
v. 5.7a (Ventana Systems). 


RESULTS 


Figure 4 shows a graphic comparison between the 
simulation results using the parameters of Table 4 
and the real data of the Hong Kong outbreak. We 
compared only the variables of the model for which 


records were found in the public databases, which 
are the same six variables shown in Figure 1. 

The simulation output for the variable ‘sick per 
day’ fit the data reported by the Hong Kong auth- 
orities (Fig. 4a), suggesting that the model was able 
to reproduce the epidemic curve. We observed that 
the number of new cases per day obtained in the simu- 
lation grew during the first 46 days. From this time, 
the number of new infections gradually fell to values 
<1 at later stages of the outbreak. 

As a consequence, the auxiliary variable that stores 
the cumulative SARS cases showed a characteristic 
sigmoidal growth (Fig. 45), consistent with the real 
data. We observed that the number of SARS cases 
grew from one at the beginning of the epidemic 
to around 1800 at the final stage of the outbreak, simi- 
lar to the 1755 cases reported by the authorities. 
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Fig. 5. Basic reproductive number obtained from the 
simulation. Ro was calculated as the ratio between the 
contagion rate and the recovery rate. 


Moreover, we observed a high fit between the model 
predictions and the real data. 

As expected, the number of deaths and recovered 
people in the simulation also grew with a sigmoidal 
shape to reach values similar to those found in the 
public databases (Fig. 4c, d). However, we observed 
a partial fit of these stock variables to the data 
reported during the outbreak. These slight mismatches 
were also observed in the flow variables for recovered 
and daily deaths (Fig. 4e, f) that may be due to delays 
in reporting of cases by the authorities. 

In short, looking the simulation results of Figure 3 
we can conclude that our model is able to reproduce 
largely the most important indicators of the SARS 
epidemic that occurred in Hong Kong in 2003. 

Focusing on the evolution of the basic reproduction 
number (Fig. 5), we note that during the early stages 
of the epidemic, Rp >1, which is consistent with the 
observed disease spread. Moreover, after day 30, Ro 
starts to drop to values <1, probably due to the im- 
plementation of containment measures by the Hong 
Kong authorities after the issuance of the first global 
alert against SARS on 12 March 2003. These results 
are consistent with a previous report showing the 
basic reproductive numbers for different SARS epi- 
demic curves, which supports the notion that our 
model is able to largely replicate the disease outbreak 
in Hong Kong [31]. 

The results of the univariate sensitivity analysis are 
shown in Figure 6. We focused our attention on the 
epidemic wave although the analysis is possible in 
other variables, as shown in the multivariate analysis 
of Figure 7. Variations in case fatality, threshold of 
cumulative attack rate and disease duration induced 
little changes in the epidemic curve (Fig. 6—c), while 
more extensive alterations in the epidemic wave were 


observed after changes in infectivity, daily contacts 
and incubation period (Fig. 6d-/). Variations in the 
case-fatality parameter do not alter the output of the 
variable sick per day (Fig. 6a), although other vari- 
ables from the model such as recovered and dead 
are highly impacted (data not shown). Small changes 
in the shape and the maximum of the curve are 
observed after modification of the parameters ‘thresh- 
old of cumulative attack rate’ and ‘disease duration’ 
(Fig. 6b, c). By contrast, changes in infectivity, daily 
contacts and to a lesser extent in incubation period 
significantly alter both the position and height of the 
maximum of the epidemic curve (Fig. 6d—/). 

The results of the multivariate sensitivity analysis 
are shown in Figure 7, we analysed the output of 
four variables: sick per day, infected, recovered, and 
dead (Fig. 7a-d, respectively). Variations in the 
model parameters clearly change the shape of the epi- 
demic curve, altering both the position and the height 
of the maximum of this variable (Fig. 7a). Similarly, 
the output of the variables infected and recovered is 
highly impacted by the changes in parameters carried 
out in the multivariate sensitivity analysis (Fig. 7), c). 
The alterations of variable outputs are clearly exem- 
plified by the variable ‘dead’ (Fig. 7d). Certain 
changes in the parameters of the model can explain 
an increase in the total number deaths during the epi- 
demic rising from about 300 to about 700. The 
observed variations in model output when the value 
of the parameters is changed support the idea that 
this model might be able to simulate different scenar- 
ios and epidemic conditions. 


DISCUSSION 


System dynamics modelling has been successfully ap- 
plied to study complex public health issues such as the 
design of optimal policies in healthcare [32], the im- 
pact of public health intervention in different situa- 
tions [33, 34], and to study disease epidemiology 
[35-37]. In the latter, system dynamics technology 
has become a powerful tool to understand and predict 
the impact of infectious diseases. Epidemiological 
models can help health authorities to make recom- 
mendations regarding intervention to fight the spread 
of directly transmissible pathogens, especially when 
empirical data is limited. In this sense, mathematical 
models have been previously used to advise health 
policies against diseases such as pandemic influenza 
[38, 39] and SARS [31, 40, 41]. 
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Fig. 6. Sensitivity of the model to changes in parameters. Univariate analysis with 1000 simulations was performed for 
each parameter: (a) case fatality, (b) threshold of cumulative attack rate, (c) disease duration, (d) infectivity, (e) daily 
contacts, and (f) incubation period. The parameters were varied considering ranges and distributions shown in Table 5. 
The solid black line represents the simulation output and the grey area represents the 95% confidence bounds. 


Several models studying SARS transmission and 
interventions have been published. These are detailed 
hybrid stochastic and compartmental models that 
successfully explain the behaviour of the epidemic 
[14, 31, 42]. Here, trying to follow Ockham’s razor 
principle, we have built a simpler deterministic 
model, which is also able to reproduce the behaviour 
of the epidemic based on its natural history and the in- 
tervention measures taken in Hong Kong. The use of 
a less complicated model could be helpful in under- 
standing the disease epidemic and also facilitating its 
reuse under other conditions. 

The epidemiological models depend on the consist- 
ency of the chosen parameters. Therefore, the accurate 
quantification of these parameters is critical to esti- 
mate the path of a disease, to predict the impact of 
possible interventions, and to inform planning and 


decision making. Here, we have combined reported in- 
formation from the SARS epidemic with an iterative 
optimization process to set the final values for 
the model parameters. Under these conditions, the 
model output fits the epidemic curve observed in the 
Hong Kong SARS-CoV outbreak (Fig. 4). 

Of the factors that influence the dynamics of infec- 
tious diseases, the person-to-person contact pattern 
has been shown to be essential in disease spread [43]. 
Our model takes this essential factor into account 
through the auxiliary variable ‘frequency of contacts’, 
which is dependent on the auxiliary variable ‘cumulat- 
ive attack rate’ and the parameters ‘daily contacts’ 
and ‘threshold of cumulative attack rate’ (Fig. 2). A 
previous work showed that mixing patterns and con- 
tact characteristics were remarkably similar across 
the different European countries analysed in that 
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Fig. 7. Multivariate sensitivity analysis of the model. A Monte-Carlo sensitivity analysis with 1000 simulations was 
performed and the results for (a) sick per day, (b) infected, (c) recovered, and (d) dead are shown. The values of the six 
parameters were altered at the same time following the distributions described in Table 5. The solid black line represents 
the simulation output and the grey area represents the 95% confidence bounds. 


study even though the average number of contacts 
recorded differed. Interestingly, the authors suggests 
that the results may well be applicable to other coun- 
tries with similar social structures, and that the initial 
epidemic phase of an emerging infection in susceptible 
populations, such as SARS was in 2003, is likely to be 
very similar [30]. During the SARS outbreak, health 
authorities, hospitals, and the overall population pro- 
gressively implemented quarantine and protection 
policies to prevent the transmission of the disease 
[15]. With this in mind, we made the assumption that 
the intervention of the health authorities caused a de- 
crease in the frequency of contacts, which in turn led a 
decrease in the rate of contagion. To mimic this event 
in our model, the number of daily contacts is regulated 
by means of a repression Hill function, which has been 
used to simulate repressor activities in complex biologi- 
cal systems [21]. 

The basic reproductive number is a key epidemiolo- 
gical variable that characterizes the potential of a dis- 
ease to spread. Several works have estimated that 
prior to the first global alert the basic reproductive 
number for SARS was >1, correlating with an expo- 
nential increase in the number of cases. However, 
the implementation of effective control measures, 
such as quarantine, isolation, and strict hygiene prac- 
tice in hospitals led a sudden decrease in Ro [14, 31]. 


The fact that, in our model, Ro drops to values <1 
around the date of the global alert reinforces the 
idea that a decrease in the frequency of contacts is 
able to effectively simulate the effects of the control 
measures established in the first stages of the epidemic. 
It is important to note that estimations of Ro in pre- 
vious publications are considerably lower (around 
2-4) than ours, which is almost 12 at the beginning 
of the outbreak [14, 31, 42]. Nevertheless, the esti- 
mated value of Rp also differs in these works and 
the credible intervals surrounding these deterministic 
estimations were wide, reaching superior values of 
almost 8. This high variability can be explained 
in part by the superspreading events that occurred in 
SARS epidemics. Superspreading is an unusual situ- 
ation, in which a single individual directly infects a 
large number of other people that has a large influence 
on the early course of the epidemic [42]. Interestingly, 
the fact that our model does not explicitly account 
these events could partially explain the very high esti- 
mated value of Ro at the beginning of the outbreak. 
The reliability of the model parameters is supported 
by univariate and multivariate sensitivity analyses 
(Figs 6 and 7). Furthermore, sensitivity analysis is a 
powerful tool to analyse the influence of certain de- 
cision making in the evolution of the epidemic. 
Taken together, these findings strongly suggest that 


our model, together with other system dynamics mod- 
els can be used by epidemiologists to investigate the 
likely consequences of future re-emergences of 
SARS-CoV based on analysis of the previous known 
epidemics. In addition, by adapting the key para- 
meters of these models or with a little change in the 
model structures, they can be used to face emerging 
outbreaks of infectious diseases, such as the recent 
MERS-CoV epidemic. In this regard, there are several 
similarities and differences which should be taken into 
account when using this model. Both SARS-CoV and 
MERS-CoV may cause severe respiratory failure, 
extrapulmonary features such as diarrhoea and also 
mild or asymptomatic cases. In contrast with SARS, 
MERS has lower human-to-human transmission po- 
tential, affects predominantly older people with 
more comorbid illness and has a higher case-fatality 
rate [44]. These factors would affect the output of 
the model variables (e.g. epidemic curve, cumulative 
cases, number of dead, etc.). 
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